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Image Restoration with Signal-dependent Camera Noise 

Ayan Chakrabarti and Todd Zickler 



Abstract — This article describes a fast iterative algorithm for image 
denoising and deconvolution with signal-dependent observation noise. 
We use an optimization strategy based on variable splitting that adapts 
traditional Gaussian noise-based restoration algorithms to account for 
the observed image being corrupted by mixed Poisson-Gaussian noise 
and quantization errors. 

Index Terms — Image Restoration, Signal-dependent noise, Poisson 
Noise, Variable Splitting. 



I. Introduction 

IMAGE restoration refers to the recovery of a clean sharp image 
from a noisy, and potentially blurred, observation. Most state- 
of-the-art restoration techniques [ ]-[ ] assume that observed image 
is corrupted by signal-independent additive white Gaussian noise 
(AWGN), since this makes both analysis and estimation significantly 
more convenient. However, observation noise in an image captured 
by a real digital camera is typically signal-dependent, and therefore 
image restoration algorithms based on AWGN models make sub- 
optimal use of the information available in the observed image. For 
example, one source of noise in recorded intensities is the uncertainty 
in the photon arrival process, which is Poisson distributed and has 
a variance that scales linearly with the signal magnitude. Therefore, 
using an AWGN model fails to account for the noise variance at 
darker pixels being lower than that at brighter ones, and can lead to 
over- smoothing in darker regions of the restored image. 

Due to these reasons, the development of restoration algorithms 
based on accurate noise models has been an area of active re- 
search [61— [10]. However, this task is made challenging by the 
fact that state-of-the-art restoration methods use sophisticated image 
priors that are defined in terms of coefficients of some spatial 
transform, and combining these priors with a non-Gaussian noise 
model significantly adds to complexity during estimation. Denoising 
methods for signal-dependent noise are either based on a (sometimes 
approximate) statistical characterization of this noise in transform 
coefficients [ ]-[ ], or use a variance-stabilization transform that 
renders the noise Gaussian followed by traditional AWGN denoising 
techniques [ ], [ ]. Recently, Harmany et al. [11] presented an 
iterative deconvolution algorithm that approximates a Poisson-noise 
likelihood cost at each iteration with a quadratic function based on 
the curvature of the likelihood at the current estimate. This technique 
was combined with various sparsity-based priors in [ ] to enable 
restoration in the presence of Poisson noise. 

We describe an iterative image restoration framework that accounts 
for the statistics of real camera noise. It is also an iterative technique 
like the one in [ ], but uses an optimization strategy known as 
variable -splitting leading to significant benefits in computational 
cost. The use of this strategy for image restoration dates back to 
the work of Geman et al. [ ], and has been used to develop 
fast deconvolution algorithms that use the so-called total-variation 
(TV) image prior model [ ], with extensions that consider non- 
quadratic penalties on the noise — including a Poisson likelihood 
cost [13]. In this paper, we deploy this technique to enable efficient 
estimation with a likelihood cost that models observation noise as a 
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combination of Gaussian noise, signal-dependent Poisson shot-noise, 
and digitization errors. Moreover, we develop a general framework 
that is not tied to any specific choice of image prior, and this allows 
us to adapt any state-of-the-art AWGN restoration technique (such as 
BM3D [ ] for denoising) for use with the proposed signal-dependent 
noise model. We demonstrate the efficacy of this approach with 
comparisons to both AWGN and signal-dependent state-of-the-art 
restoration methods. 

II. Observation Model 

Let x(n) be the latent noise-free sharp image of a scene corre- 
sponding to an observed image y(n), where n £ R 2 indexes pixel 
location. We assume that a spatially-uniform blur acts on the scene, 
with a known kernel k(n). We let Xk(n) = (x * k)(n) denote the 
blurred image that would have been observed in the absence of noise. 
The observation y(n) can then be modeled as 



y(n) = Q [y(n)\ , y(n) = y k (n) + z(n), 



(1) 



where Q(-) is a quantization function used to digitize the analog 
sensor measurement y(n), which we in turn model as the sum of a 
scaled Poisson random variable yk(n) with mean Xk(n), and zero- 
mean Gaussian random noise z(n) with variance a 2 . We examine 
the statistical properties of each component of the model in (1), and 
define a likelihood function based on these properties with the aim 
of enabling accurate yet efficient inference of x(n) from y(n). 

A. Shot Noise 

At each location n, the random variable yk captures the uncertainty 
in the photon arrival process, and is modeled with a scale parameter 

a > as ayk ~ V(axk), i.e., 



P(yk\x k ) 



(ax k ) c 



(2) 



(otykY- 

The difference between the observed y k and its mean Xk is referred 
to as shot noise, and has a signal-dependent variance equal to Xk/ot. 
The parameter a depends on the ISO setting of the camera, and a 
high value of a corresponds to a low ISO setting and a higher signal- 
to-noise ratio (SNR). We define Lp(xk;yk, ol) as the corresponding 
negative log-likelihood (up to a constant) of the observed pixel yk 
being due to Xk'. 



L P (x k ;y k ,a) = ax k - ay k log x k . 



(3) 



Note that Lp is a convex function of Xk (since d 2 Lp/dx 2 z — 
ayk/xl > 0, Vxfc), with a minimum at Xk = yk- 

B. Gaussian Noise 

In addition to shot noise, the measurement y may be corrupted by 
other signal-independent noise sources, such as thermal and amplifier 
noise, which we model with the AWGN variable z ~ A/"(0, a 2 ) [ ], 
[10]. Combining this with the Poisson model in (2), we have: 



p(y\x k ) oc ^2 

r=0 



(ax k ) r e 



- exp 



(y - r/a) 2 
2a 2 



(4) 
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Unfortunately, the above expression can not be computed in closed 
form, and therefore we employ an approximation to define the log- 
likelihood function for this case. We note that y is a mixed Poisson- 
Gaussian random variable with mean Xk and variance (xfc/a + a 2 ), 
and approximate it with a shifted Poisson likelihood as 



Lpc(xk;y, a, o) = L P (xk + aa 2 ;y + aa 2 ;a). 



(5) 



C. Quantization 

Finally, we account for the observed intensity y being quantized by 
a function Q(-) that maps every interval y+] to a single value y. 
We consider the general case where quantization is possibly preceded 
by a non-linear map for gamma-correction: 



y - 



(6) 



where g corresponds to the gamma-correction exponent (1 for linear 
data, and typical 2.2 for sRGB), and |_'_U denotes rounding off 
intervals of width q. Note that y here denotes a linearized version of 
the camera image, i.e., one where the inverse of the gamma-correction 
function has been applied to the quantized observations. The interval 
may therefore be asymmetric around y and have variable 
widths, and the mean and variance of y given y are given by 

(y^o + q /2Y +1 -(y 1 ^ -q/2Y +1 
m q (y) = - V J V ' 



q(9 + 1) 



y 1/9 + 



a/2) 



2g+l 



1/9 



a/2) 



2g+l 



q (2 9+ l) -<(V).0) 

Note that this variance is only signal-dependent when g / 1, since 
for g = 1, Cq(y) — q 2 /12 is independent of y. We incorporate these 
to obtain the overall likelihood function as 



L(x k ;y,a,cr, Q) = Lp G (x k ;m q (y),a, yj 'a 2 + cr 2 (y)). 



(8) 



III. MAP-based Restoration 

With a suitably defined prior model p(x) for natural images, we 
recover the maximum a-posteriori (MAP) estimate x(n) of x(n) from 
y(n) as: 



x — argmin<3>(x) + L((x * k)(n)]y(n)), 



(9) 



where = — logp(x), and L(-) is the likelihood function defined 
in (8) (with the arguments a, a, Q omitted for brevity). 

The main obstacle to computing a solution for (9), even in the 
absence of blur (i.e., k = 5), is the fact that natural image priors 
are best defined in terms of coefficients in some transform domain. 
Unlike the AWGN case where noise in the coefficients of the 
observed image is also independent and Gaussian, estimation under 
the noise model in (8) is challenging because of the complexity in 
characterizing the statistics of noise in the transform domain. 

A. Variable Splitting 

We use an optimization approach similar to the one in [13], that 
allows us to deal with minimizing the prior and likelihood terms in 
their respective natural domains, i.e., in terms of transform coeffi- 
cients and individual pixels respectively. Specifically, we recast the 
unconstrained minimization problem in (9) as an equality-constrained 
optimization with the addition of new variables r(n): 

x = arg min 4>(x) + L(r(n); y(n)), 



subject to: r(n) = (x * k)(n). 



(10) 



Clearly, the problems in (9) and (10) are equivalent. However, this 
modified formulation can be solved using an efficient iterative ap- 
proach that allows the noise in each pixel to be treated independently. 

Note that instead of choosing a specific image prior, we assume 
the existence of a baseline AWGN-based image restoration algorithm 
that (perhaps implicitly) defines and also provides an estimator 

function G(-): 



G(y, k, a) = arg min 4>(x) + 



[y{n) - (x * k){n)f 
2a 2 



(11) 



While we treat G(-) as a "black box" in general, the appendix 
describes a "parallel" optimization algorithm for a special case when 
the baseline algorithm is itself based on variable- splitting. 

B. Minimization with Augmented Lagrangian Multipliers 

We use the augmented Lagrangian multiplier method [13] to solve 
this constrained optimization, and define a new augmented cost 
function that incorporates the equality constraint: 



C(x,t,X) = $(x) + 



P 



53(r(n) -x k {n)) 2 



■^A(n) (r(n) - x k (n)) 



+ ^L(r(n);y(n)), (12) 



where Xk(n) = (x * k)(n), A(n) are the Lagrange multipliers, and 
P is a fixed scalar parameter. Note that the Lagrangian terms in 
this expression are augmented with an additional quadratic cost on 
the equality constraint. This additional cost is shown to speed up 
convergence, and since it has a derivative equal to zero when the 
equality constraint is satisfied, it does not affect the Karush-Kuhn- 
Tucker (KKT) condition at the solution. 

The solution to (10) can be reached iteratively by making the 
following updates to x(n),r(n) and A(n) at each iteration: 



T 



arg min C (x, r* , A* ), 

X 

arg min C(x t+1 , r, A*), 



(13) 
(14) 



\\n) - 7 (n)/?(r t+1 (n) - x\ + \n)), (15) 



where x t , r* , 7* refers to the values of these variables at iteration t, 
and the step size j(n) lies between and j max = (y/E + l)/2. 

The update to x in (13) involves minimizing the sum of the prior 
term <3>(x) with a uniformly- weighted quadratic cost, and this can be 
achieved using the baseline estimator G(-) as 



= G(r\n)-p- 1 X\n),k,p- 1 ) 



(16) 



The update to r in (14) involves a minimization that is independent 
of the prior term, and can be done on a per-pixel basis. For each n, 
we solve for dC/dr(n) = to obtain 



r t+i 



(n) = (6+ V6 2 + 4c)/2, 



<v [a 2 + a 2 (y(n))) , 



6= (^(rO+ZT^n)) - a/5- 1 , „ qK 

c = a (a 2 + ^(2/(n))) (x* +1 (n) + /T^ra)) + ar'^feW). 

(17) 

C. Algorithm Details 

We begin the iterations with r(n) = m q {y{n)) and A(n) = 0, 
and stop when the relative change in x falls below a threshold: 

||x t+1 (n) -x*(n)" 2 



< e. 



(18) 
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We find that it is optimal to vary the step-size 7(72) at each pixel 
(but this is kept fixed across iterations), with a higher step-size 
for pixels with higher observed intensities y(n) that have a higher 
expected noise variance. Specifically, we vary the step-size linearly 
with respect to y(n), between 7 max /2 and 7 max . 

Finally, the choice of /3 involves a trade off between accuracy 
and speed of convergence. We find that it is best to choose a value 
inversely proportional to the average expected noise variance in the 
image, i.e. P = /?o/o"avg» where 

aavg = a -1 avg{ra q (j/(n))} + a 2 + avg{ag (y(n))}, (19) 

and /3o depends on the choice of the baseline algorithm G(-). 

IV. Experimental Results 

We use synthetically blurred and noisy images to compare the 
proposed approach to traditional AWGN-based methods for de- 
convolution, as well as Poisson-Gaussian noise-based methods for 
denoising [ ], [ ]. Table I shows deconvolution performance on 
three standard images blurred with circular pill-box kernels (that 
correspond to typical defocus blur) of different radii r, in the presence 
of Poisson-Gaussian noise with different values of a and a, as well as 
with quantization errors (marked as +Q in the table) corresponding 
to 8-bit quantization post gamma-correction (q = 1/256, g = 2.2). 

We report performance in terms of PSNR (averaged over ten in- 
stantiations of the noise) for the AWGN-based ADM-TV method [ ], 
and for the proposed method using the ADM-TV method as the 
baseline. We use the parallel optimization approach as described in 
the appendix, and use the same baseline parameters (Pv = 200, k = 
8,e=10 -5 )in both cases. We set f3o = 16 for the proposed method, 
and use the mean noise- variance cr| V g for the AWGN results. We 
note that the proposed approach offers a distinct advantage over the 
baseline method with better estimates in all cases. Figure 1 shows 
an example of the observed and restored images, and we see that the 
proposed method is able to account for the noise variance at darker 
pixels being lower, yielding sharper estimates in those regions. Our 
approach typically requires only two to three times as many iterations 
as the baseline method to converge. 

Next, we show results in Fig. 2 for the two demonstration examples 
provided by authors of [ ], each blurred with a 5 x 5 kernel and 
corrupted only by Poisson noise. We show the running times and 
PSNR values of the restored images, from both the proposed method 
(using the same parameters as above), and the SPIRAL technique [11] 
with a TV prior (which yields the highest PSNR values amongst the 
different priors). In both cases, our approach yields estimates with 
higher accuracy, and offers a significant advantage in computational 
cost — for the cameraman example, our method converges in just 
25 iterations, while the SPIRAL technique requires a 100 iterations, 
with each iteration in-turn calling an iterative baseline TV- solver. 

Finally, we report denoising performance (i.e., k — S) in Table. II, 
for the standard test cases used in [10] with Poisson-Gaussian noise. 
The state-of-the-art AWGN techniques for denoising tend to be more 
sophisticated than those for deconvolution. We use the BM3D [ ] 
algorithm, which uses a complex adaptive image prior, as the baseline 
estimator for our approach in this case (with Pq = 2,e = 10" 3 ). In 
addition to PSNR values for the proposed method, we show results 
for the baseline AWGN method (with cr| V g as the noise variance), and 
for the Poisson-Gaussian denoising algorithms described in [ 0] and 
[8]. Note that the method in [ ] also uses BM3D as a baseline. We 
use the same notation as in [ ] to describe the noise parameters in 
Table II, where noise is synthetically added by first scaling the input 
image to a certain peak value, instantiating Poisson random variables 
with these scaled intensities, and then adding Gaussian noise with 




Input (PSNR=18.64dB) 




AWGN Method [ ] (21.48 dB) Prop. Method (22.29 dB) 



Fig. 1. Deconvolution results on the Camerman image, blurred with a circular 
kernel of radius 9 pixels with noise parameters a = 1024, a = 10 -4 , q = 
1/256, and g = 2.2. Note that the proposed method correctly accounts for 
the noise variance being lower at darker pixels, and recovers sharper estimates 
in those regions in comparison to the baseline AWGN method of [4]. 

variance a 2 . The reported PSNR values are again averaged over ten 
instantiations of the noise. Figure 3 shows an example of input and 
restored images for this case. 

We find that the proposed approach is competitive with these 
existing methods, with the highest PSNR in a majority of the test 
cases. However, it is important to remember that our approach is 
iterative, while the algorithms in [8], [10] are single-shot. Table. II 
reports the mean number of iterations required in each case, and 
we see that convergence is usually quick, requiring roughly three to 
seven calls to the baseline method. We also show results for denoising 
using the simple TV prior (which was used for the deconvolution 
results in Table. I), and note that using the BM3D method as the 
baseline instead leads to a significant improvement. This highlights 
the importance of the flexibility that our approach offers, in allowing 
the use of any baseline AWGN restoration method rather than being 
tied to a fixed prior. 

V. Conclusion 

In this paper, we introduce a framework for image restoration in 
the presence of signal-dependent noise. We describe an observation 
model that accounts for camera sensor measurements being corrupted 
by both Gaussian and signal-dependent Poisson noise sources, and 
for errors from the subsequent digitization of these measurements. 
Then, we use variable-splitting to derive a fast iterative scheme that is 
able to adapt existing AWGN-based restoration methods for inference 
with this noise model. A MATLAB implementation of the algorithm, 
along with scripts to generate the results presented in this paper, is 
available for download at http://vision.seas.harvard.edu/PGQrestore/. 

The flexibility in being able to incorporate any AWGN-based 
method as a baseline means that we can draw on a considerable 
amount of existing research on image statistics. Our approach can 
therefore be easily used for restoration of any class of images (such 
as medical images, or astronomical data) which is corrupted by 
noise with similar characteristics, by using an appropriate class- 
specific AWGN-restoration technique as the baseline. Moreover, the 
optimization scheme described here can be adapted to other noise 
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TABLE I 

Deconvolution performance (PSNR in dB) for Poisson-Gaussian and Quantization 
Noise (indicated with +Q), and blur with circular pill-box kernels of radius r 



cm 


(7 


Method 


Cameraman 


Lena 


Boats 


r = 5 


r = 7 


r = 9 


r = 11 


r = 5 


r = 7 


r = 9 


r = 11 


r = 5 


r = 7 


r = 9 


r = 11 


1024 


io- 4 


Input 


20.55 


19.46 


18.64 


17.97 


25.02 


23.62 


22.58 


21.76 


23.12 


22.05 


21.35 


20.81 


AWGN [ ] 


23.21 


22.27 


21.50 


20.97 


28.41 


26.99 


26.23 


25.40 


25.55 


24.29 


23.64 


23.00 


Proposed 


24.08 


23.09 


22.29 


21.69 


28.70 


27.34 


26.51 


25.70 


26.03 


24.80 


24.12 


23.42 


1024 


io- 1 


Input 


17.23 


16.69 


16.24 


15.85 


18.54 


18.19 


17.87 


17.57 


18.27 


17.89 


17.61 


17.38 


AWGN [ ] 


21.05 


20.18 


19.41 


18.73 


25.57 


24.53 


23.67 


22.95 


23.16 


22.42 


21.86 


21.41 


Proposed 


21.55 


20.60 


19.82 


19.11 


26.11 


24.96 


24.06 


23.22 


23.57 


22.73 


22.10 


21.57 


1024 


io- 4 
+Q 


Input 


20.55 


19.46 


18.64 


17.97 


25.02 


23.61 


22.58 


21.76 


23.11 


22.05 


21.34 


20.81 


AWGN [ ] 


23.21 


22.26 


21.50 


20.96 


28.41 


26.99 


26.23 


25.39 


25.55 


24.29 


23.64 


23.00 


Proposed 


24.08 


23.08 


22.28 


21.68 


28.70 


27.34 


26.50 


25.70 


26.02 


24.79 


24.11 


23.41 


1024 


io- 1 
+Q 


Input 


17.38 


16.81 


16.35 


15.94 


18.58 


18.21 


17.89 


17.58 


18.33 


17.93 


17.64 


17.41 


AWGN [ ] 


20.97 


20.13 


19.38 


18.70 


25.56 


24.53 


23.66 


22.94 


23.14 


22.40 


21.85 


21.40 


Proposed 


21.50 


20.56 


19.80 


19.09 


26.10 


24.96 


24.06 


23.22 


23.56 


22.72 


22.09 


21.56 


256 


io- 4 


Input 


19.93 


18.97 


18.23 


17.62 


23.29 


22.30 


21.51 


20.86 


21.96 


21.13 


20.54 


20.10 


AWGN [ ] 


22.20 


21.26 


20.50 


19.87 


27.26 


25.90 


25.05 


24.18 


24.48 


23.40 


22.73 


22.18 


Proposed 


22.88 


21.93 


21.12 


20.51 


27.64 


26.28 


25.50 


24.60 


24.95 


23.76 


23.05 


22.43 


256 


io- 1 


Input 


16.93 


16.42 


15.99 


15.62 


18.09 


17.76 


17.47 


17.20 


17.86 


17.52 


17.26 


17.04 


AWGN [ ] 


20.96 


20.09 


19.32 


18.65 


25.42 


24.39 


23.55 


22.84 


23.03 


22.33 


21.79 


21.34 


Proposed 


21.47 


20.52 


19.75 


19.04 


25.99 


24.85 


23.96 


23.12 


23.46 


22.64 


22.01 


21.51 


256 


io- 4 
+Q 


Input 


19.93 


18.97 


18.23 


17.62 


23.29 


22.29 


21.51 


20.85 


21.96 


21.12 


20.54 


20.09 


AWGN [ ] 


22.19 


21.26 


20.50 


19.86 


27.26 


25.89 


25.05 


24.18 


24.47 


23.40 


22.73 


22.18 


Proposed 


22.88 


21.93 


21.12 


20.51 


27.64 


26.28 


25.50 


24.60 


24.95 


23.76 


23.05 


22.43 


256 


io- 1 
+Q 


Input 


17.08 


16.54 


16.10 


15.71 


18.14 


17.80 


17.50 


17.23 


17.92 


17.56 


17.30 


17.07 


AWGN [ ] 


20.88 


20.04 


19.29 


18.62 


25.40 


24.38 


23.55 


22.84 


23.01 


22.31 


21.77 


21.33 


Proposed 


21.41 


20.47 


19.72 


19.01 


25.98 


24.84 


23.95 


23.12 


23.45 


22.63 


22.00 


21.50 


64 


io- 4 


Input 


18.06 


17.43 


16.90 


16.45 


19.65 


19.19 


18.79 


18.42 


19.09 


18.64 


18.31 


18.03 


AWGN [ ] 


21.35 


20.45 


19.68 


18.99 


25.95 


24.83 


23.96 


23.18 


23.41 


22.61 


22.02 


21.56 


Proposed 


21.98 


20.99 


20.20 


19.58 


26.51 


25.31 


24.44 


23.56 


23.87 


22.95 


22.31 


21.77 


64 


io- 1 


Input 


15.88 


15.49 


15.15 


14.84 


16.64 


16.40 


16.19 


15.98 


16.51 


16.26 


16.06 


15.90 


AWGN [ ] 


20.63 


19.81 


19.04 


18.43 


24.88 


23.95 


23.16 


22.54 


22.64 


22.01 


21.53 


21.11 


Proposed 


21.19 


20.29 


19.50 


18.79 


25.55 


24.48 


23.58 


22.86 


23.11 


22.36 


21.79 


21.33 


16 


io- 4 


Input 


14.26 


13.99 


13.77 


13.55 


14.50 


14.35 


14.21 


14.08 


14.42 


14.25 


14.13 


14.02 


AWGN [ ] 


20.10 


19.31 


18.62 


18.15 


24.02 


23.20 


22.55 


21.99 


21.99 


21.50 


21.08 


20.67 


Proposed 


20.87 


19.97 


19.22 


18.54 


24.90 


23.90 


23.05 


22.40 


22.58 


21.93 


21.45 


21.01 


16 


io- 1 


Input 


13.22 


13.01 


12.83 


12.66 


13.33 


13.23 


13.12 


13.02 


13.36 


13.23 


13.13 


13.04 


AWGN [ ] 


19.73 


18.97 


18.41 


18.03 


23.48 


22.76 


22.18 


21.65 


21.64 


21.20 


20.81 


20.40 


Proposed 


20.45 


19.60 


18.84 


18.26 


24.35 


23.42 


22.67 


22.06 


22.20 


21.63 


21.20 


20.76 


4 


io- 4 


Input 


9.02 


8.97 


8.92 


8.87 


8.72 


8.69 


8.65 


8.60 


8.81 


8.77 


8.74 


8.70 


AWGN [ ] 


18.31 


17.99 


17.71 


17.41 


20.86 


20.45 


20.07 


19.67 


19.85 


19.57 


19.26 


18.99 


Proposed 


18.86 


18.38 


18.08 


17.74 


21.70 


21.15 


20.72 


20.25 


20.53 


20.18 


19.79 


19.43 


4 


io- 1 


Input 


8.69 


8.64 


8.60 


8.55 


8.39 


8.35 


8.32 


8.28 


8.49 


8.45 


8.42 


8.39 


AWGN [ ] 


18.23 


17.93 


17.65 


17.36 


20.64 


20.24 


19.87 


19.49 


19.70 


19.43 


19.14 


18.89 


Proposed 


18.62 


18.22 


17.92 


17.58 


21.34 


20.83 


20.41 


19.96 


20.26 


19.94 


19.56 


19.24 



models, as long as the likelihood functions for those models are 
convex, or can be so approximated. 

Appendix 

The ADM-TV method [ ] is a popular choice for deconvolution 
under AWGN, and is itself based on variable- splitting. When using 
this method as the baseline algorithm, it is possible to adopt an 
approach where its internal optimization is effectively done in parallel 



to that for the noise model in our framework. We describe this 
approach in detail in this appendix. The ADM-TV method uses an 
image prior that penalizes the TV-norm of image gradients: 

Q(x) = k £ fe|(x*V0(n)| 2 , (20) 

n y i 

where V* are gradient niters, and k is a model parameter. The 
MAP estimate in this case is computed by introducing additional 
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Fig. 2. Deconvolution performance with Poisson Noise in comparison to the SPIRAL-TV method [ ]. The running times are indicated in brackets. 

TABLE II 

Denoising performance (PSNR in dB) with Poisson-Gaussian Noise 



Image 


Peak 




!N"oisy 






PTTRF-T FT TR1 


JT 1 UU . TLJ 1V1 J VJ 


7T ±lCIClllUll& 


Prnn +TV 

-T 1 UU . I 1 V 


Cameraman 


1 


0.1 


3.20 


18.50 


20.23 


20.35 


20.71 


7.2 


17.42 


2 


0.2 


6.12 


20.95 


21.93 


21.60 


22.12 


5.1 


18.53 


5 


0.5 


9.83 


23.55 


24.09 


23.33 


24.10 


4.0 


21.08 


10 


1 


12.45 


25.10 


25.52 


24.68 


25.57 


4.0 


22.88 


20 


2 


14.76 


26.50 


26.77 


25.92 


26.81 


4.0 


24.55 


30 


3 


15.91 


27.10 


27.30 


26.51 


27.30 


4.0 


25.26 


60 


6 


17.49 


27.97 


28.07 


27.35 


28.10 


3.0 


26.12 


120 


12 


18.57 


28.52 


28.57 


27.89 


28.61 


3.0 


26.56 


Fluorescent 
Cells 


1 


0.1 


7.22 


19.64 


24.54 


25.13 


24.91 


8.8 


20.79 


2 


0.2 


9.99 


22.24 


25.87 


26.25 


26.16 


7.3 


22.68 


5 


0.5 


13.37 


25.43 


27.45 


27.60 


27.72 


5.5 


25.21 


10 


1 


15.53 


27.53 


28.63 


28.59 


28.77 


5.0 


27.11 


20 


2 


17.21 


29.18 


29.65 


29.47 


29.69 


5.0 


28.37 


30 


3 


17.97 


29.91 


30.16 


29.84 


30.18 


4.0 


28.87 


60 


6 


18.86 


30.72 


30.77 


30.42 


30.75 


4.0 


29.38 


120 


12 


19.39 


31.15 


31.14 


30.70 


31.09 


4.0 


29.60 


Lena 


1 


0.1 


2.87 


19.87 


22.59 


22.83 


22.58 


7.0 


17.01 


2 


0.2 


5.82 


22.58 


24.34 


24.16 


24.14 


5.1 


19.68 


5 


0.5 


9.54 


25.42 


26.17 


25.74 


26.19 


4.0 


23.18 


10 


1 


12.19 


27.20 


27.72 


27.27 


27.80 


4.0 


25.41 


20 


2 


14.53 


28.66 


29.01 


28.46 


29.12 


4.0 


27.05 


30 


3 


15.72 


29.42 


29.69 


29.12 


29.75 


3.5 


27.69 


60 


6 


17.35 


30.37 


30.51 


29.91 


30.57 


3.0 


28.12 


120 


12 


18.48 


30.98 


31.05 


30.51 


31.11 


3.0 


28.08 



auxiliary variables di(n) corresponding to the image gradients, and 
the estimation problem is cast as: 



The iterative optimization for this case proceeds as follows: 



x = argmm^^ \di(n)\ 2 + L(r(n); y(n)), 

n V i 

subject to: di(n) = (x * Vi)(ra); r(n) = (x * k)(n). 



in [13]: 

C(x, di, r, A) = k J^2\Mn)\ 2 + 

n V i 

- ^2 ^i( n )(di(n) - Xi(n)) + 

n,i 





d? 1 * 


- arg min C(x t , di, r*. A*), 


(23) 




x t+1 <• 


- argminC(x, c/- +1 , r*, A*), 

X 


(24) 


(21) 


T t+1 * 


. J^f t+1 »t+l x t\ 
- arg mm C(x ,a { , r, A ) , 

T 


(25) 


case as 


A t+1 (n) * 


- A t (n)- 7 (n)/3(r f+1 (n)-4 +1 (n)), 


(26) 




A^ +1 (n) <■ 


- A*(n) - %™/3v{dl +1 {n) -xl +1 (n)). 


(27) 



2 

2 



y^(di(n) - Xi(n)) 2 
^(r(n) -x k (n)f 



+ L(r(n);y(n)), (22) 



where Xi(n) = (x * Vi)(n), and the second and third term above 
encode the gradient equality constraint in (21). 



Note that this is essentially equivalent to the optimization framework 
in (13)-(15), with (23), (24) corresponding to the x update step in 
(13), and (27) being an extra update step for the additional Lagrange 
multipliers (using a fixed step size equal to 7 max , as in [ ]). The 
solution to the r update step is identical to the one described in (17) 
in Sec. Ill, since the minimization in (25) for this update depends 
only on the last three terms in the cost in (22). 

The update steps to d\ and x in (23), (24) are largely independent of 
the noise model, and are similar to those described in [ ]. Specifically, 
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Fig. 3. Denoising results (with inset close-ups) on the Fluorescent cells 
image, with Poisson-Gaussian noise corresponding to a peak input intensity 
of 5 and a = 0.5. 



each d\ +1 (n) can be computed independently on a per-pixel basis: 



jt+i 



(n) 



d(n) 



d{n) = ^ \xj(n) +/3 v 1 A-(r 



(28) 



The updated value of x can then be computed efficiently in the Fourier 
domain as: 



gv EjZ[V£jMW] +;8J r [<!]V[f(n) 
/3v EJ^[V l ]| 2 +/3|^[fc]| 2 



(29) 



where, 

4(n) = " ^AJW, f(n) = T*(r») - /T^n), (30) 



and J 7 and J 7-1 refer to the forward and inverse discrete Fourier 
transforms. To account for the periodicity assumption with the Fourier 
transform in (29), we extend di and f in each direction by six times 
the size of the kernel k before computing the Fourier transform, with 
zeros for di, and values that linearly blend the intensities at opposite 
boundaries for f. After computing and we crop 

them back to their original sizes. 

We begin the iterations with r(n) = m q (y(n)), X(n) = 
0, x(n) = and Xi(n) = 0, and choose the gradient filters V* to 
correspond to horizontal and vertical finite-differences (i.e., [—1, 1]). 
We vary the step sizes 7(71) as described in Sec. III. Moreover, we 
do not apply the updates in (25) and (27) to r(n) and A(n) for the 
first few iterations (six in our implementation), during which time the 
algorithm proceeds identically to the AWGN ADM-TV [ ] method 
with input m q (y(n)) and noise variance 
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